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^ ■ Abstract: We present a phenomenological approach (EPOS), based on the par- 



ton model, but going much beyond, and try to understand proton-proton and 



deuteron-gold collisions, in particular the transverse momentum results from all 
the four RHIC experiments. It turns out that elastic and inelastic parton ladder 
splitting is the key issue. Elastic splitting is in fact related to screening and sat- 
uration, but much more important is the inelastic contribution, being crucial to 
understand the data. We investigate in detail the rapidity dependence of nuclear 
effects, which is actually relatively weak in the model, in perfect agreement with 
the data, if the latter ones are interpreted correctly. 

1 Introduction 

Interesting new results have been observed in heavy ion collisions at RHIC, a 
large fraction of which being based on transverse momentum spectra. High 
transverse momenta seem to be suppressed [1,2], systematically different for 
different hadron species [3] . 

But any quantitative expression of a suppression or an enhancement needs a 
reference, and here one usually refers to proton-proton scattering. Rather than 
working with the spectrum, one investigates the ratio Raa of the nucleus-nucleus 
[AA] spectrum to the proton-proton [pp) result, the so-called nuclear modification 
factor, 

R _ 1 dn AA dnP' p ^ 
A iV co ii d 2 p t dy d?p t dy ' 
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The normalization factor l/-/V coll has been chosen such that at large transverse 
momenta one expects R AA to become unity. Here, dn refers to the number of 
produced particles per inelastic interaction. 

The first problem is therefore to understand sufficiently well proton-proton 
scattering. This is far from trivial. Experimentally, it is difficult to really ac- 
cess the full inelastic cross section, the interaction triggers tend to miss a more 
or less large fraction of the events. Theoretically, proton-proton is far from be- 
ing fully understood, apart from perturbative calculations concerning very large 
transverse momenta. 

A second problem arises due to the fact that even being sure about the ob- 
servation of a non-trivial behavior of Raa, we want to know whether this effect 
is really a collective one, providing evidence of the formation of a quark gluon 
plasma, and not something we observe already in proton-nucleus. This was the 
main purpose to study, in addition to gold-gold [Au Au), as well deuteron-gold 
[dAu) collisions at RHIC [4-8], with quite interesting results: the strong high p t 
suppression in AuAu seems to be absent in dAu, so we have clearly a final state 
effect. 

Many features of dAu seem to be qualitatively understood, employing the sat- 
uration model [9-12], a recombination model [13], an improved parton model 
[14-16], the AMPT model [17]. But, what is really missing is a global and quan- 
titative investigation: can we understand ALL the data presented so far by ALL 
the experiments, for pp and dAu, in a single approach. This gives also the op- 
portunity to cross-check the different experiments, which is not so obvious to do 
directly. 

The purpose of this paper is to present a phenomenological approach (EPOS), 
based on the parton model, but going much beyond, and try to understand pp 
and dAu, as far as the transverse momentum results from all the four RHIC 
experiments are concerned. It turns out that parton splitting (or better parton 
ladder splitting) is the key issue, which is related to screening and saturation, 
but there are other important consequences, which are crucial to understand 
the data. 

2 Improved Parton Model with Remnants 

The new approach we are going to present is called EPOS, which stands for 

D Energy conserving quantum mechanical multiple scattering approach, 
based on 

D Partons (parton ladders) 

d Off- shell remnants 

D Splitting of parton ladders 

We are going to explain the different items in the following (the parton splitting 
will be discussed in a later section). 
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Figure 1 : Elementary parton-parton scattering: the hard scattering in the middle 
is preceded by parton emissions (initial state space-like cascade); these partons 
being usually off-shell, they emit further partons (final state time-like cascade). 
For all this we use a symbolic parton ladder. 

One may consider the simple parton model to be the basis of hadron-hadron 
interaction models at high energies. It is well known that the inclusive cross 
section is given as a convolution of two parton distribution functions with an 
elementary parton-parton interaction cross section. The latter one is obtained 
from perturbative QCD, the parton distributions are deduced from deep inelastic 
scattering. Although these distributions are taken as black boxes, one should 
not forget that they represent a dynamical process, namely the successive emis- 
sion of partons (initial state space-like cascade), which have to be considered in a 
complete picture. In addition, the produced partons are generally off-shell, giving 
rise again to parton emissions (final state time-like cascade). All this is sketched 
in fig. [H where we also indicate that we refer to this whole structure as "parton 
ladder", with a corresponding simple symbol, to simplify further discussion. 

For practical calculations, each parton ladder is finally translated into two 
color strings, which fragment into hadrons. This is a purely phenomenological 
procedure for the non-perturbative hadronization process. 

Actually our "parton ladder" is meant to contain two parts: the hard one, 
as discussed above, and a soft one, which is a purely phenomenological object, 
parametrized in Regge pole fashion, for details see Appendix iBl 

Still the picture is not complete, since so far we just considered two interact- 
ing partons, one from the projectile and one from the target. These partons leave 
behind a projectile and target remnant, colored, so it is more complicated than 
simply projectile/target deceleration. One may simply consider the remnants to 
be diquarks, providing a string end, but this simple picture seems to be excluded 
from strange antibaryon results at the SPS [18]. 

We therefore adopt the following picture, as indicated in fig. El not only 
a quark, but a two-fold object takes directly part in the interaction, being a 
quark- antiquark, or a quark-diquark, leaving behind a colorless remnant, which 
is, however, in general excited (off-shell). So we have finally three "objects", 
all being white: the two off-shell remnants, and the parton ladders between 
the two active "partons" on either side (by "parton" we mean quark, antiquark, 
diquark, or antidiquark) . We also refer to "inner contributions" (from parton 
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Figure 2: The complete picture, including remnants. The remnants are an im- 
portant source of particle production at RHIC energies. 




Figure 3: Inner contributions, from the parton ladder (full lines), and "outer" 
contributions, from the remnants (dashed lines), to the rapidity distribution of 
hadrons. (Artists view) 



ladders) and "outer contributions" (from remnants), reflecting the fact that the 
remnants produce particles mainly at large rapidities and the parton ladders at 
central rapidities, see fig. [3J. Whereas the outer contributions are essentially 
energy independent, apart of a shift in rapidity, the inner contributions grows 
with energy, to eventually dominate completely central rapidities. But at RHIC, 
there is still a substantial remnant contribution at mid-rapidity. 
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Figure 4: The two elements of the multiple scattering theory: open ladders, 
representing inelastic interactions, and closed ladders, representing elastic in- 
teractions. 

Even inclusive measurements require often more information than just in- 
clusive cross sections, for example via trigger conditions. Anyhow, for detailed 
comparisons we need an event generator, which obviously requires information 
about exclusive cross sections (the widely used pQCD generators are not event 
generators in this sense, they are generators of inclusive spectra, and a Monte 
Carlo event is not a physical event). This problem is known since many years, the 
solution is Gribov's multiple scattering theory, employed since by many authors. 
This formulation is equivalent to using the eikonal formula to obtain exclusive 
cross sections from the knowledge of the inclusive one. 

We indicated recently inconsistencies in this approach, proposing an "energy 
conserving multiple scattering treatment". The main idea is simple: in case 
of multiple scattering, when it comes to calculating partial cross sections for 
double, triple ... scattering, one has to explicitly care about the fact that the 
total energy has to be shared among the individual elementary interactions. 

A consistent quantum mechanical formulation requires not only the consid- 
eration of the (open) parton ladders, discussed so far, but also of closed ladders, 
representing elastic scattering, see fig. |4j The closed ladders do not contribute 
to particle production, but they are crucial since they affect substantially the 
calculations of partial cross sections. Actually, the closed ladders simply lead to 
large numbers of interfering contributions for the same final state, all of which 
have to be summed up to obtain the corresponding partial cross sections. For 
details see appendix 

We can do the complicated calculations, since we fit for example the result of 
a numerical calculation of a squared amplitude corresponding to a (open) parton 
ladder of energy y^i, using a simple form as 13 , which allows then to perform an- 
alytical calculations. Furthermore, we employ very sophisticated Markov chain 
techniques to generate configurations according to multidimensional probability 
distributions. For details see appendix |0 an d appendix iDl 

Important concerning numerical results: There are a couple of parameters 
which determine the parameterization of the soft elementary interaction (soft 
Pomeron), which are essentially fixed to get the pp cross sections right. The 
pQCD parameters (soft virtuality cutoff, K-factor, parton emission cutoff, parton- 
hadron coupling) are fixed to provide a reasonable parton distribution function 
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Figure 5: Basic parton-parton interaction in nucleus-nucleus collisions: a pro- 
jectile parton always interacts with exactly one parton on the other side, elasti- 
cally (closed parton ladder) or inelastically (open parton ladder) . 

(which we calculate, it is not input!). 

We assume the remnants to be off-shell with probability p Q , a mass distribu- 
tion given as 

proboc ikr 2a °, (2) 

with parameter values which are not necessarily the same for diffractive and non- 
diffractive interactions (the latter ones being defined to be those without parton 
ladders). We use currently for p Q 0.75 (dif) and 0.95 (nondif), and for a 0.75 (dif) 
and 1.1 (nondif). Those excitation exponents may give rise to quite high mass 
remnants, RHIC and also SPS data seem to support this. High mass remnants 
will be treated as strings 

There are four important fragmentation parameters: the break probability 
(per unit space-time area) p B , which determines whether a string breaks earlier 
or later, the diquark break probability p D , the strange break probability p s , and 
the mean transverse momentum p t of a break, with obvious consequences for 
baryon and strangeness production, and the p t of the produced hadrons. We use 
three sets of these parameters, for the three types of strings: soft-, kinky(hard)-, 
remnant- strings. We do not really use the full freedom of these parameters, but 
one single set would not work - if we are interested in high precision. Somewhat 
surprising: p s is 0.14 for soft and 0.06 for kinky strings. Maybe this reflects the 
fact that soft strings may have low masses, where strangeness is suppressed, 
and which needs some compensation. The parameter p D is as well bigger for soft 
compared to kinky strings. 



3 Splitting of Parton Ladders 

Let us first consider very asymmetric nucleus-nucleus collisions, like proton- 
nucleus or deuteron-nucleus. The formalism developed earlier for pp can be 
generalized to these nuclear collisions, as long as one assumes that a projectile 
parton always interacts with exactly one parton on the other side, elastically or 
inelastically (realized via closed or open parton ladders), see fig. [5J 
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Figure 6: Inelastic and elastic "rescattering" of a parton from the parton ladder 
with a second target parton. We talk about (inelastic and elastic) splitting of a 
parton ladder. 

We employ the same techniques as already developed in the previous sec- 
tion. The calculations are complicated and require sophisticated numerical 
techniques, but they can be done. The corresponding results for dAu will be 
discussed later. 

In case of protons (or deuterons) colliding with heavy nuclei (like gold) , there 
is a complication, which has to be taken into account: suppose an inelastic 
interaction involving an open parton ladder, between a projectile and some target 
parton. The fact that these two partons interact implies that they are close in 
impact parameter (transverse coordinate). Since we have a heavy target, there 
are many target partons available, and among those there is a big chance to find 
one which is as well close in impact parameter to the two interacting partons. In 
this case it may be quite probable that a parton from the ladder interacts with 
this second target parton, inelastically or elastically, as shown in fig. |BJ 

Let us first discuss the effects of elastic splitting. The squared amplitude for 
an elementary inelastic interaction involving two partons with light cone momen- 
tum shares x + = 2p + /^ and x~ = 2p~/ v / s can be parametrized quite accurately 
as 

a (x + f(x-f, (3) 

with two parameters a and (3 depending on the squared energy s and the impact 
parameter b [*/s is the proton-proton cms energy). Any addition of an elastic 
contribution (closed ladder), be it in parallel or via splitting, provides an inter- 
ference term, contributing negatively to (partial) cross sections. So an additional 
elastic leg, even though it does not affect particle production, it provides screen- 
ing. Model calculations show that adding elastic splittings to the basic diagrams, 
modifies the corresponding squared amplitude as 

a {x + f(x-f +e , 

and therefore the whole effect can be summarized by a simple positive exponent 
e, which suppresses small light cone momenta. So the existence of many target 
partons effectively screens small x contributions, which agrees qualitatively with 
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Figure 7: Transport of transverse momentum via an attached closed ladder, 
which may be even enhanced in case of diquarks. 



Figure 8: Hadron production in case of inelastic ladder splitting. 

the concept of saturation. But this is only a part of the whole story, several other 
aspects have to be considered. 

One effect is the transport of transverse momentum via an attached closed 
ladder, as shown in fig. Such a transport we use already in the basic par- 
ton model, when it comes to diffractive scattering, realized via a closed ladder. 
Here, some transverse momentum transfer is needed to explain the transverse 
momentum spectra of protons at large x (in the diffractive region). In case of 
diffractive target excitation, the projectile gets simply a p t kick. We should have 
the same phenomenon in case of elastic splitting: the ladder parton involved in 
the interaction should get a p t kick in the same way as the proton in diffractive 
scattering. This could be even more effective for diquarks (two attached ladders), 
leading finally to baryon production. 

Let us turn to inelastic splitting, fig. |5] . Consider the example shown in the 
figure. In the upper part, there is only an ordinary parton ladder, so we expect 
"normal" hadronization. In the lower part, we have two ladders in parallel, which 
are in addition close in space, since they have a common upper end, and the 
lower ends are partons close in impact parameter, so the hadronization of the 
two ladders is certainly not independent, we expect some kind of a "collective" 
hadronization of two interacting ladders. Here, we only considered the most 
simple situation, one may also imagine three or more close ladders, hadronizing 
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collectively. 

So far we discussed in a qualitative fashion the consequences of elastic and 
inelastic parton ladder splitting. The strength of the effects will certainly depend 
on the target mass, via the number Z of partons available for additional legs. The 
number Z of available partons will also increase with energy, so at high enough 
energy the above-mentioned effects can already happen in pp collisions. 

4 Realization of Ladder Splitting Effects 

The basic quantity for a numerical treatment of the ladder splitting effects is the 
number Z of partons available for additional legs, more precisely we have a Z T for 
counting legs on the target side, and Z P for counting legs on the projectile side. 
Let us treat Z T (corresponding discussion for Z P ). Consider a parton in projectile 
nucleon i which interacts with a parton in target nucleon j. The number Z T (i,j) 
of addition legs has two contributions, one counting the legs attached to the 
same nucleon j, and one counting the legs attached to the other nucleons j' ^ j. 
We assume the following form: 

Z T (i,j) = z exp{-b 2 j /2b 2 ) + 

+ £ 4 exp(-^,/26 2 ), 

target nucleons j'^j 

where 6^ is the distance in impact parameter between i and j . The coefficients z 
and z' Q depend logarithmically on the energy, as 

z = w z In s/ s M , 

z'o = w z a/(1ii s/s m ) 2 + W}/ , 

and the impact parameter width is b = wb v^imW^. with parameters wb, wz, 
w m , and s M - We then define 

z t(j) = /^2 z T{hj)- 

i 

We suppose that all the effects of the parton ladder splitting can be treated 
effectively, meaning that the correct explicit treatment of splittings is equiva- 
lent to the simplified treatment without splittings, but with certain parameters 
modified, expressed in terms of Z. 

This is not only to simplify our life. Even an explicit dynamical treatment will 
stay a phenomenological approach, with many uncertainties, for example about 
the splitting vertices, and much more. So we prefer to have simple parameteri- 
zations rather than a very complicated but uncertain dynamical treatment. 

So which quantities depend on Z, and how? In the following the symbols a { 
are constants, used as fit parameters. The elastic splitting leads to screening, 
which is expressed by the screening exponents e = es (for soft ladders) and e = eh 
(for hard ladders), and here we assume 

e s = asP s Z, (4) 
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e H = a H p H Z , (5) 

where (3 S and (3 H are the usual exponents describing soft and hard amplitudes 
(see appendix P. 

A second effect is transport of transverse momentum, here we suppose 

Ap t = a T p n q Z, (6) 

where n q is the number of quarks of the objects in the hadronization process (1 
for quarks, 2 for diquarks), and p = 0.5 GeV is just used to define a scale. 

Let us come to the collective hadronization. We will actually "absorb" the 
multiple ladders into the remnants, which are usually treated as strings. Now 
we treat them as strings with modified string break parameters, to account for 
the collective hadronization. We modify the break probability (per unit space- 
time area) p B , which determines whether a string breaks earlier or later, the 
diquark break probability p D , the strange break probability p s , and the mean 
transverse momentum p t of a break, as 



Pb -> Pb - a B Z , (7) 

Pd -> Pd (1 + &d Z) , (8) 

Ps -> p s (l + a s Z), (9) 

p t -> p t (l + apZ), (10) 



with positive parameters a*. So with increasing Z, a reduced p B will lead to more 
particle production, an increased p D , p s , p t , will lead to more baryon-antibaryon 
production, more strangeness production, and an increased p t for each string 
break. 

The parameters s M , w it and a, are chosen to reproduce mainly RHIC pp and 
dAu data, but also the energy dependence of cross sections and multiplicities 
from SPS to Tevatron. The best fit parameters are shown in table [TJ 



5 Results for Proton-Proton 

Ladder splitting is quite important for pp at very high energies, where cross sec- 
tions and multiplicities are considerably suppressed, due to screening. At RHIC 
energies, however, the effects are small: the total cross section is reduced by 
5%, the multiplicity by 10%. Concerning the transverse momentum spectra to 
be discussed in detail in the following, the effect is hardly visible. 

In order to compare to the charged particle p t spectra in pp from the different 
experiments (STAR, PHENIX, BRAHMS), one has first to understand what has 
been measured. One wants to measure the inelastic differential yield, 

d 3 n incl 1 dV nel 
dyd 2 p t (Tine! dyd 2 pt 

where a w 42 mb is the inelastic pp cross section, and <fV ncl / 'dyd 2 p t represents 
the inclusive differential cross section for inelastic events. In practice there is an 
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Transverse momentum tran snort 


0.025 


a B 


Break parameter 


0.070 


a D 


Diquark break probability 


0. 110 


a s 


Strange break probability 


0.140 


ap 


Average break transverse momentum 


0.150 



Table 1: Best fit values for splitting parameters. We included in the fit as well 
data not shown in this paper. 

event trigger like the beam beam counter (BBC) , which only counts a fraction of 
the events, missing in particular low multiplicity events. 

UA5 [20] actually used a similar trigger to define non single diffractive (NSD) 
events. The NSD differential yield is given as 



d 3 n NSD 1 d 3 a 



NSD 

(12) 



dyd 2 p t cxnsd dyd 2 p t ' 

where a 35 mb is the NSD pp cross section, and d 3 cx NSD /dyd 2 p t represents the 
inclusive differential cross section for NSD events. To be clear: NSD is not an ab- 
solute definition, it is defined via the acceptance of the UA5 detector! In figure |9l 
we show the corresponding pseudorapidity distribution for NSD events, slightly 
higher than the one for inelastic events. For the simulation of NSD events, we use 
simply the same requirement as used in the experiment (coincidence of charged 
particles in a forward and a backward pseudorapidity interval) . 

In case of STAR, one could as well define NSD as the events accepted by the 
BBC. Let us do so for the moment, and use the term NSD BBC . What is actually 
done is somewhat different. The differential cross section is multiplied by 30/26, 
in order to correspond to what Pythia defines to be non single diffractive, cor- 
responding to 30mb. Let us call this NSD PYT . Actually the inelastic differential 
yield does not change, however, it is interpreted as spectrum for NSD PYT events. 
Then again based on Pythia, it is argued that the inelastic differential yield for 
inelastic events is obtained essentially (with a small correction at small p t ) by 
multiplying with 30/42 (just the ratio of the cross sections), since SD events do 
not contribute to particle production. So after all, the originally measured dif- 
ferential yield (referring to NSD BBC ) and the inelastic one differ essentially by a 
factor 42/30 = 1.4 - what is not at all what we observe, simulating NSD events 
with the BBC trigger condition, and comparing with inelastic events. As seen in 
fig. 1 1 01 the ratio of the NSD BBC differential yield to the inelastic differential yield, 
rather than being 1.4, differs considerably as a function of p t , and in addition 
depends on the particle species. 
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Figure 9: Pseudorapidity distribution for inelastic and NSD events in pp colli - 
sions.The lines are EPOS results, the points data [20]. The dotted lined repre- 
sents the "inner contribution" to the inelastic distribution (many particles are 
coming from remnants!). 
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Figure 10: Ratio of the NSD BBC differential yield to the inelastic differential yield, 
in pp collisions, for a pions [n), kaons (K), and protons (p). 



Based on the above discussion, we will simulate NSD BBC differential spectra 
(we actually get exactly 26 mb for the NSD BBC cross section), and compare with 
STAR'S published NSD results, since they are identical to the NSD BBC differential 
spectra. In this way we avoid all these problems related to Pythia correction 
procedures. In the following, NSD refers always to NSD BBC . 

In fig. II 1L we show p t spectra for NSD events, compared to STAR data [21], 
and inelastic events, compared to PHENIX data [6,22]. Simulation and data 
agree within 15% (over 6 orders of magnitude). This good agreement is only 
possible after our re-interpretation of NSD, see the discussion above. 

When studying (later) dAu collisions, there will be plenty of discussion con- 
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Figure 11: Differential yields in pp collisions as a function of p t , for (from top to 
bottom) charged particles (over 2) for NSD events, charged particles (over 2) for 
inelastic events, and neutral pions for inelastic events. Lines are EPOS simula- 
tions, points are data from STAR [21] and PHENIX [6,22]. The two agree within 
1 5% (over 6 orders of magnitude) . 
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Figure 12: Inelastic differential yields in pp collisions as a function of p t , for (from 
top to bottom) charged particles (over 2) at rj = 0, at rj — 1, negative particles at 
7] = 2.2, at?7 = 3.2 (always displaced by factors of 10). Lines are EPOS simulations, 
poins are data [7]. We also plot (dashed) the simulation curve at 77 = 0, multiplied 
by 0.1, 0.01, 0.001, to have a reference. 

cerning the (pseudo)rapidity dependence of certain effects. It is therefore neces- 
sary to first check the (pseudo)rapidity dependence of p t spectra for pp. In fig. 
1121 we plot inelastic differential yields as a function of p t , at different pseudora- 
pidities; 37 = 0, 77 = 1, 77 = 2.2, and 77 = 3.2 We show EPOS simulations compared to 
BRAHMS data [7]. We also plot (dashed) the simulation curve at 77 = 0, multiplied 
by 0. 1, 0.01, 0.001, to have a reference for the results at the other pseudorapidi- 
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Figure 13: The target Z as a function of centrality, expressed in terms of the 
number of binary collisions, for dAu. 
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Figure 14: The Z distribution, for different centrality classes, 
ties. The spectra clearly get softer with increasing 77. 

6 Results Deuteron-Gold 

All screening effect are linear in Z, so it is worthwhile to first investigate Z. In 
very asymmetric collisions as dAu, the projectile Z is essentially zero, whereas 
the target Z differers considerably from zero. As show in fig. 1131 (and obvious 
from the definition) Z T increases linearly with the number of collisions. So Z is 
essentially a centrality measure. In fig. 1141 we show the Z distribution for the 
different centrality classes. In this way one understands easily how the different 
centrality classed are affected by the splitting effects. 

Although we are mainly interested here in transverse momentum spectra, 
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Figure 15: Pseudorapidity spectra in minimum bias dAu collisions. Lines are 
EPOS simulations, points are data from PHOBOS [23] (circles), STAR [4] (tri- 
angles), BRAHMS [24] (squares). We also show the inner and the outer target 
contributions to the simulated distribution. 

we still show first of all the pseudorapidity spectra, which finally determine the 
normalization of the p t spectra. In fig. 1151 we show pseudorapidity spectra in 
minimum bias dAu collisions: EPOS simulations, compared to data from from 
PHOBOS [23], STAR [4], and BRAHMS [24]. We also show different contribu- 
tions to the simulated distribution. We distinguish inner and outer (projectile 
and target) contributions, where the outer contributions are meant to contain 
the multiple ladders, originating from ladder splittings, treated in a "collective" 
way, as discussed above. The inner contribution comes from ordinary ladders in 
the middle. The asymmetry of the distribution is clearly due to to the target con- 
tribution (the projectile contribution, not shown, is very small). Also not shown 
here is the result for "no splitting"; for central collisions, ladder splitting leads 
to an overall reduction of dn/drj of about 30%. This is due to the fact that first 
of all the two splitting effects "screening" and "string break delay" are relatively 
small, and secondly work in opposite direction. In figs. 1161 and 1171 we show 
pseudorapidity spectra for central and peripheral dAu collisions. 

Let us now turn to p t spectra. One of the first observations concerning p t 
spectra in dAu collisions was the fact that not only the nuclear modification 
factor shows a non-trivial behavior, but this behavior seems also to be strongly 
pseudorapidity dependent, even when varying r/ only by one unit. We want to 
investigate this question in the following. 

In fig. 1181 we show transverse momentum spectra of charged particles in 
dAu collisions at different centralities and at different pseudorapidities. The four 
figures represent minimum bias, central (0 — 20%), mid-central (20 — 40%) , and 
peripheral (40 — 100%) collisions. For each figure, spectra for four pseudorapidity 
intervals are shown: [—1,-0.5], [—0.5,0], [0,0.5], [0.5,1]. We simply refer to the 
corresponding mean values, r] = —0.75, i] = —0.25, r] = 0.25, r] = 0.75. For better 
visibility, the different curves have been displaced by factors of 10. Solid lines 
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Figure 16: Pseudorapidity spectra for central dAu collisions. Solid lines are EPOS 
simulations for 0-20%, dashed lines are simulations for 0-30%. Points are data 
from PHOBOS [25] (circles), STAR [4] (triangles), BRAHMS [24] (squares). 
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Figure 17: Pseudorapidity spectra for peripheral dAu collisions. Solid lines 
are EPOS simulations for 40-100%, dashed lines are simulations for 60-80%. 
Points are data from PHOBOS [25] (circles), STAR [4] (triangles), BRAHMS [24] 
(squares). 
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Figure 18: Transverse momentum spectra of charged particles in dAu collisions 
at different centralities and at different pseudorapidities. The four figures rep- 
resent minimum bias, central (0 - 20%), mid-central (20 - 40%) , and peripheral 
(40 - 100%) collisions. For each figure, from top to bottom: rj = -0.75, rj = -0.25, 
7] = 0.25, i] = 0.75. Solid lines are EPOS simulations, points are data [4]. The 
different curves have been displaced by factors of 10. 



are EPOS simulations, points are data [4], both agree within 10-20%. Although 
looking directly at spectra does not really allow to see systematic differences 
between the different curves, it is still useful to first check that the absolute 
curves agree, before investigating ratios. 

To observe any "anomalous" behavior, one usually plots ratios, like the nu- 
clear modification factor, defined earlier. The disadvantage is the fact that the 
corresponding pp spectrum has to be known, with a sufficient precision. An 
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Figure 19: R cp ratios at different pseudorapidities (77 = —0.75, 77 = —0.25, 77 = 0.25, 
77 = 0.75). Solid lines are EPOS simulations, points are data [4]. The dotted lines 
are EPOS simulations, with parton ladder splitting turned off. 

alternative procedure is the use of ratios of central to peripheral results, 

1 J„ central i A^. peripheral 

° P N coT al d 2 Ptdy 1 yyp^pheraj d 2 p t dy 

In fig. 1191 we show the R cp ratios at different pseudorapidities (77 = —0.75, 
77 = —0.25, 77 = O.2577 = 0.75). Here, central refers to — 20% and peripheral to 
40 — 100%. Solid lines are EPOS simulations, points are data [4]. We also show 
the corresponding EPOS results, with parton ladder splitting turned off. These 
curves are cut off at p t = 3 GeV/c, to avoid that the strong statistical fluctua- 
tions spoil the figure. The no-parton-ladder -splitting curve increases slowly with 
p t , in the shown range it stays well below one. It will finally reach one. The 
full EPOS simulations show quite a different behavior, the ratio R cp increases 
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Figure 20: Ratio of charged particle spectra at backward to forward pseudora- 
pidities {q = —0.75/ r] = 0.75), in minimum bias dAu collisions. Solid lines are 
EPOS simulations, points are data [4]. The dotted lines are EPOS simulations, 
with parton ladder splitting turned off. 

strongly between 1 and 2 GeV/c, to stay constant (or decrease) beyond. The sta- 
tistical fluctuations do not really allow very precise predictions beyond 4 GeV/c. 
The strong increase between 1 and 2 GeV/c is due to collective hadronization on 
the target side, which leads to an increased transverse momentum production. 
As can be seen from fig. 1151 target side hadronization extends even to forward 
pseudorapidities, so it is quite visible in the whole r/ range [-1,1]. The effect 
is simply somewhat stronger at backward compared to forward rapidity, since 
target hadronization contributes more. But the difference is not so big. The in- 
crease of R cp with p t is partly also due to the momentum transfer from the elastic 
splitting, which should affect equally backward and forward pseudorapidities. 
The variation of the shape of R cp with pseudorapidity is quite small, the main 
modification is actually an overall factor due to the fact, that the particle density 
increases towards smaller pseudorapidities, as seen from fig. 1151 

A direct way to investigate the pseudorapidity dependence of spectra is pro- 
vided by the ratio of spectra at backward to forward pseudorapidities, like 77 = 
—0.75/ 77 = 0.75, as shown in fig. 1201 Here, one observes a slight increase between 
and 2 GeV/c. This means that the R cp at backward pseudorapidity increases 
a bit more than the one at forward pseudorapidity, which we understand such 
that there is somewhat more target side collective hadronization at backward 
pseudorapidity. 

We now consider an even larger pseudorapidity variation: we investigate how 
the nuclear modification factors vary in the pseudorapidity range to 3.2. Before 
comparing to data, we show the results of full EPOS simulations, as well as those 
with parton ladder splitting turned off. In fig. 1211 we show the nuclear modifi- 
cation factors for charged particles in minimum bias dAu collisions, at different 
pseudorapidities: 77 = 0, 77 = 1, 77 = 2.2, 77 = 3.2. Whereas the no-splitting curves 
hardly change with p t , and decrease with pseudorapidity, the full calculations 
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Figure 21: Nuclear modification factors RdAu for charged particles in minimum 
bias dAu collisions, at different pseudorapidities: rj = 0, 77 = 1, 77 = 2.2, 77 = 3.2 
(from top to bottom). The solid lines are full EPOS simulations, the dotted lines 
are EPOS with parton ladder splitting turned off. 

show of course the same decrease with pseudorapidity, but all curves increase 
substantially with p t between 1 and 3 GeV/c. This confirms the observation 
already made earlier by studying the variation in the 77 range - 1 to 1 . 

In the following we will compare the simulations with data from all the four 
RHIC experiments. In fig. 1221 we collect all published data on charged particle 
nuclear modification factors in minimum bias dAu collisions at (or close to) 77 = 0, 
together with the corresponding simulations. We show minimum bias results at 
7] = from STAR [5], at 77 = 0.4 from PHOBOS [8], 0-88% centrality results at 
77 = from PHENIX [6], and minimum bias data at 77 = from BRAHMS [7]. We 
also show minimum bias EPOS simulations at 77 = 0, at 77 = 0.4, not feed down 
corrected minimum bias simulations, and 0-88% centrality results at 77 = 0. 
We first of all observe that the different simulation results are quite close to 
each other, so changing slightly the pseudorapidity, the centrality definition, 
doing or not feed down corrections, does not affect the final result too much. 
The variation of the experimental data is much bigger. On the upper end we 
have the STAR data, but based on the above discussion on pp results, we expect 
that the pp reference spectrum is 10-20% too low, which means RdAu is 10-20% 
too high. The corresponding reduction would bring the STAR data down to the 
EPOS simulation curve (full line), and agree with the PHENIX data, and with 
the PHOBOS data. BRAHMS is on the lower end, but within the error bars 
compatible with the simulation curve. 

In fig. 1231 we consider charged particle nuclear modification factors in min- 
imum bias dAu collisions at (or close to) 77 = 1 , together with the corresponding 
simulations. We show minimum bias data at 77 = 0.8 from PHOBOS [8], at 77 = 1 
from BRAHMS [7]. We also show minimum bias EPOS simulations at 7/ = 1, at 
77 = 0.8, both being very close to each other. The data are somewhat lower, but 
the curves are within the error bars. A systematic difference may again be due 
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Figure 22: Nuclear modification factors RdAu for charged particles in minimum 
bias dAu collisions at (or close to) rj = 0. The different lines are full EPOS simu- 
lations: minimum bias at 77 = (full), at 77 = 0.4 (dashed-dotted), not feed down 
corrected (dotted), 0-88% centrality, at 77 = (dashed). The points are minimum 
bias data at rj = from STAR [5] (rhombs), at 77 = 0.4 from PHOBOS [8] (squares), 
0-88% centrality data at 77 = from PHENIX [6] (circles), minimum bias results 
at 77 = from BRAHMS [7] (triangles). 

to the pp reference. To investigate this, we also plot a "mixed" RdAu'- the nuclear 
spectrum is taken from BRAHMS, but we use the EPOS pp reference. The result 
(squares) now exceeds the simulation curve. 

In fig. 1241 we finally compare EPOS simulations and data from BRAHMS [7] 
at 77 = 2.2 and 77 = 3.2. Data and simulations agree quite well. 

7 Summary 

In this paper, we have presented a phenomenological approach, called EPOS, 
based on the parton model, but going much beyond. There are two very im- 
portant "nuclear effects": elastic and inelastic parton ladder splitting, which in 
principle occurs already in pp scattering, but which becomes really visible when 
systems with large numbers of partons (like nuclei) are involved. 

Elastic splitting is in fact related to screening and saturation, but much more 
important is the inelastic contribution, being crucial to understand the data. 

The main effect (at least concerning the observables investigated in this paper) 
is due to the fact that inelastic splitting (bifurcation of parton ladders) , leads to a 
modified hadronization process, a "collective hadronization" of multiple, parallel 
parton ladders, on the target side, in case of dAu. This is the equivalent of string 
fusion, if one uses the language of strings. But contrary to the usual string 
fusion picture, here we do not have complete ladders which behave collectively, 
but only the bifurcated ones on the target side. 
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Figure 23: Nuclear modification factors -R<l4« for charged particles in minimum 
bias dAu collisions at (or close to) rj = 1. The different lines are full EPOS simula- 
tions: minimum bias at rj — 1 (full), at rj = 0.8 (dashed). The points are minimum 
bias data at rj = 0.8 from PHOBOS [8] (rhombs), at rj = 1 from BRAHMS [7] (cir- 
cles) . The squares represent RdAu calculated from the BRAHMS dAu p t spectrum 
and the pp simulation result. 




Figure 24: Nuclear modification factors R dAu for charged particles in minimum 
bias dAu collisions at r] = 2.2 (left) and r] = 3.2 (right). The lines are full EPOS 
simulations, the points are data from BRAHMS [7]. 
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Concerning p t spectra, the main effect of the collective hadronization is a 
p t broadening. This is certainly what is needed, but real evidence for our pic- 
ture can only come from a very detailed comparison with all corresponding data 
currently available. For this purpose we considered all published nuclear modi- 
fication factor data, concerning charged particles, from all the four RHIC experi- 
ments. 

We investigated in detail the rapidity dependence of nuclear effects, which is 
actually relatively weak in the model, in perfect agreement with the data, if the 
latter ones are interpreted correctly, and if one considers really ALL available 
data, and not just a convenient subset. 
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A Appendix on Multiple Scattering Theory 



We first consider inelastic proton-proton scattering. We imagine an arbitrary 
number of elementary interactions to happen in parallel, where an interaction 
may be elastic or inelastic. The inelastic amplitude is the sum of all such contri- 
butions, 




with at least one inelastic elementary interaction involved. To calculate cross 
sections, we need to square the amplitude, which leads to many interference 
terms, like for example 




which represents interference between the first and the second diagram of the 
above-mentioned sum of terms. Using the above notations, we may represent 
the left part of the diagram as a cut diagram, conveniently plotted as a dashed 
line: 



The squared amplitude is now the sum over many such terms represented by 
solid and dashed lines. 

When squaring an amplitude being a sum of many terms, not all of the terms 
interfere - only those which correspond to the same final state. For example, a 
single inelastic interaction does not interfere with a double inelastic interaction, 
whereas all the contributions with exactly on inelastic interaction interfere. So 
considering a squared amplitude, one may group terms together representing 
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Figure 25: Class of terms corresponding to one inelastic interaction. 




Figure 26: Class of terms corresponding to two inelastic interactions. 

the same final state. In our pictorial language, this means that all diagrams with 
one dashed line, representing the same final state, may be considered to form 
a class, characterized by m — 1 - one dashed line ( one elementary cut) - and 
the light cone momenta x + and x~~ attached to the dashed line (denning energy 
and momentum of the Pomeron). In fig. 1251 we show several diagrams belonging 
to this class, in fig. 1261 we show the diagrams belonging to the class of two 
inelastic interactions, characterized by m = 2 and four light-cone momenta xf, 
Generalizing these considerations, we may group all contributions 
with m inelastic interactions [m dashed lines = m cut Pomerons) into a class 
characterized by the variable 

We then sum all the terms in a class K, 

Q(K) = ^{all terms in class K}. 
The cross section is then simply a sum over classes, 




Q depends implicitly on the energy squared s and the impact parameter b: Q = 
Q,( s,b \ The individual terms j d 2 bVt(K), represent partial cross sections, since they 
represent distinct final states. They are referred to as topological cross sections. 

The above concepts are easily generalized to nucleus-nucleus scattering, an 
example for a diagram representing a contribution to the squared amplitude is 
shown in fig. 12 71 We may also define classes, which correspond to well defined 
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final states, in our notation a given number of dashed lines between nucleon 
pairs. We may number the pairs as 1, 2, 3, ... k ... , AB. We define m k to be 
the number of inelastic interactions (cut Pomerons) of the pair number k. The 
fi th of these m k cut Pomerons is characterized by light cone momenta xt, x7. So 
a class may be characterized by 

K = {m k ,x^,x^}. 

We sum all terms in a class to obtain again a quantity called f^ s,6 )(lT), such that 
the cross section can be written as a sum over classes 



as in the case of proton-proton scattering. Here, however, b is a multidimensional 
variable representing the impact parameter b and the transverse distances b k of 
all the nucleon-nucleon pairs, 

b = {b ,b 1 ,...,b AB }, (14) 
and / d 2 b is a symbolic notation, meaning 

A B 

fd 2 b= [ d\ I f[d%?T A (b?)f[d 2 bfT B (bf), (15) 

J ^ J i=l j=l 

with A, B being the nuclear mass numbers and with the so-called nuclear thick- 
ness function T A (b) being defined as the integral over the nuclear density pa(b), 

T A (b) := J dz PA (VWT^). (16) 
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One can prove 

K 

which is a very important result justifying our interpretation of Q(K) to be a 
probability distribution for the configurations K. This provides also the basis for 
applying Monte Carlo techniques. 

The function Q is the basis of all applications of this formalism. It provides the 
basis for calculating (topological) cross sections, but also for particle production, 
thus providing a consistent formalism for all aspects of a nuclear collision. 

For the sake of completeness, we provide the formula for Q, expressed in 
terms of the elementary interactions G fc (indices k express that G depends on the 
considered pair because of ladder splitting correction) and some vertex function 
F [26]: 



AB ( 1 m k -\ 

X+ X ) = U\—^U s ' b ^ ^ (x*«*,x**, s, b) , (17) 

fc=l I™*- M =l J 



with a function $ representing summations over all uncut contributions corre- 
sponding to a given cut configuration, given as 



$^(^,^,5,6) = E-E (18) 

AB f l k ^ AB ( l k 

n \ n d3t t\ d3t k,x \ n \ n n _g *( 5 m> ^ ^ & * 

k=l Ia=1 J fc=l I fc ' A=l 



TT f 

I 1 remn 
i=l 

The different arguments are denned as 




X+ 


~ { X k,fj} > 


x- 


= i X k,n) ' 


x wo) 


= {xT J }: 


^targ 


= {^r}' 


m 


= {™fc}, 



and 



targ 



(19) 
(20) 

(21) 

(22) 
(23) 



= I" E <„ t 24 ) 

7r(fe)=i 

= 1 - E • (25) 

r(fc)=j 



The functions 7r(/c) and r(fc) refer to the projectile and the target nucleons par- 
ticipating in the k th interaction (pair k) and F remn is the vertex function to the 
projectile or target remnant. 
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B Appendix on Soft and Hard Elementary Interac- 
tions 



What are actually these G functions, which represent elementary interactions? 
To explain this we have to discuss some basic facts about scattering. Let T be 
the elastic scattering amplitude T for a 2 — > 2 scattering. The 4-momenta p and 
p' are the ones for the incoming particles , p = p + q and p' = p' — q the ones for 
the outgoing particles, and q the 4-momentum transfer in the process. We define 
as usual the Mandelstam variables s and t. Using the optical theorem, we may 
write the total cross section as 

<7tot(s) = ^2ImT(M = 0). (26) 
We define the Fourier transform f of T as 

f(s,b) = d 2 q L e-^ b T(s,t), (27) 

using t = —q\, and we define G as 

G(s,b) = j-2lmf(s,b). (28) 

One can easily verify that 

<7tat(s) = Jd 2 bG(s,b), (29) 

which allows an interpretation of G(s, b) to be the probability of an interaction at 
impact parameter b. 

This are actually such "G functions" which appear in the multiple scatter- 
ing formulas. The fact that an elementary scattering of nucleon constituents 
(quarks) is composed of soft and hard components (and a semihard one, as dis- 
cussed later) can be expressed as 

G = G SO ft + C hard . (30) 



Soft Interactions 

Let us first consider a soft interaction (with only small virtualities involved) . We 
assume that nucleon constituents (quarks, diquarks,...) from projectile and tar- 
get may interact with the result of producing many (low p t ) hadrons. We refer to 
this as an elementary soft interaction of nucleon constituents. The correspond- 
ing amplitude (T-matrix element) is assumed to be of the form 

T so{t (s, t) = i 87r s 7^ art i — J exp(A soft (s/ s ) t) , (31) 
\ so{ {z) = 2R 2 + a' soft \nz, (32) 
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with parameters a soft , a' soft , 7 part , R 2 part , and a scale s = lGeV 2 . We define G so{t (s, b) 
to be twice the imaginary part of the Fourier transform (with respect to k± = -t) 
of T soft , as discussed above. We find 

G soit (s, b) = : \ (-) QS ° ft 1 ex P f- Tx ^ / • (33) 



Koi{s/so) \s J \ 4A soft (s/s ) 



Hard interactions 

We now proceed to the case of hard scattering , being the other extreme, when 
all internal intermediate partons are characterized by large virtualities Q 2 > Ql. 
Here, the corresponding amplitude T^ ard (s, t) of the scattering of two partons with 
flavors j and k can be calculated using the perturbative QCD techniques [27, 
28]. In the leading logarithmic approximation of QCD, summing up terms where 
each (small) running QCD coupling constant a s (Q 2 ) appears together with a large 
logarithm 1ii(<5 2 /Aq CD ) (with \q C d being the infrared QCD scale), and making use 
of the factorization hypothesis, one obtains the contribution of the corresponding 
cut diagram for t = q 2 = as the cut parton ladder cross section a^ rd (s, Ql) 
(strictly speaking, one obtains the ladder representation for the process only 
using axial gauge), where all horizontal rungs of the ladder are the final (on- 
shell) partons and the virtualities of the virtual t-channel partons increase from 
the ends of the ladder towards the largest momentum transfer parton-parton 
process. We have 

<Us,Ql) = ^2ImT^ rd M = 0) 

/A n ml 
dx + B dx- B dpl-^{x + B x B s )P l) (34) 

x E^ CD (x B ,Ql,M 2 )E^B,QlM 2 )e(M 2 -Q 2 ) , 

Here da^ l OTn / dp\ is the differential 2^2 parton scattering cross section, p\ is 
the parton transverse momentum in the hard process, m, I and x% are corre- 
spondingly the types and the shares of the light cone momenta of the partons 
participating in the hard process, and Mp is the factorization scale for the pro- 
cess (we use Mp = p\/4). The 'evolution function' E^ D (Ql, Mp, z) represents the 
evolution of a parton cascade from the scale Ql to Mp, i.e. it gives the number 
density of partons of type m with the momentum share z at the virtuality scale 
Mp, resulted from the evolution of the initial parton j, taken at the virtuality scale 
Ql. The evolution function satisfies the usual DGLAP equation [29-32] with the 
initial condition E^ D (Ql,Ql, z) = 5 J m 5(1 — z). The factor K ~ 1.5 takes effectively 
into account higher order QCD corrections. 

In the following we shall need to know the contribution of the uncut parton 
ladder T^ d (s, t) with some momentum transfer q along the ladder (with t = q 2 ). 
The behavior of the corresponding amplitudes was studied in [33] in the leading 
logarithmic(l/x ) approximation of QCD. The precise form of the corresponding 
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amplitude is not important for our application; we just use some of the results 
of [33], namely that one can neglect the real part of this amplitude and that it 
is nearly independent on t, i.e. that the slope of the hard interaction Rl ard is 
negligible small, i.e. compared to the soft Pomeron slope one has Rl ard — 0. So 
we parameterize T^. d (s,t) in the region of small t as [34] 

TlUs, t) = is at d (s, Ql) exp (R^ d t) , (35) 

The corresponding "G function" is obtained by calculating the Fourier transform 
f of T and dividing by the initial parton flux 2s, 

GiUs,b) = y s 2lmf£ rd ^b), (36) 

which gives 

G hL( s > & ) = J d2 Q± exp(-iq±b) 2ImT h J a fc rd (s, -q]_), 

or explicitly 

G - (s - 6) -^ (s - gS )«b exp (-i)' 1371 

This is not yet the hard G function we are looking for, there is still one ele- 
ment missing, namely the relation between the partons j and k and the incident 
nucleons. There are two elements relevant: the vertex function, and the soft 
pre -evolution. 



The vertex function 

The coupling of the multiple (say n) elementary interactions to the nucleons is 
expressed via projectile and target nucleon vertex functions, having the form 



1 - £ x k exp -Rl IK-t^) > (38) 

\ k=l / \ k=l J k=l 

with 

F£M = In z~ Qpart , (39) 



= x a — e(x)G(l-x), (40) 

with parameters R%, j N , a part , a^ mn . All parameters with an "N" sub (super) script 
are specific for projectile/target nucleons and may be different for other kinds 
of hadrons, like pions, kaons. The arguments x k and ql ± are respectively the 
light cone momentum fractions and squared transverse momenta of the nucleon 
constituents participating in an interaction. All factors apart of F remn can be 
absorbed into the elementary G functions. 
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Soft pre-evolution 

In case of sea quarks and gluons being at the end of a parton ladder, the 
momentum share x\ of the "first" parton is typically very small, leading to an 
object with a large mass of the order Q\jx\ between the parton and the pro- 
ton [35]. Microscopically, such 'slow' partons with x\ <C 1 appear as a result of 
a long non-perturbative parton cascade, where each individual parton branch- 
ing is characterized by a small momentum transfer squared Q 2 < Q\ and nearly 
equal partition of the parent parton light cone momentum [36,37]. When cal- 
culating proton structure functions or high -p t jet production cross sections, this 
non-perturbative contribution is usually included into parameterized initial par- 
ton momentum distributions at Q 2 = Qq. However, the description of inelastic 
hadronic interactions requires to treat it explicitly in order to account for sec- 
ondary particles produced during such non-perturbative parton pre-evolution, 
and to describe correctly energy-momentum sharing between multiple elemen- 
tary scatterings. As the underlying dynamics appears to be identical to the 
one of soft parton-parton scattering considered above, we treat this soft pre- 
evolution as the usual soft emission. We account for this by introducing "soft 
pre -evolutions functions" E so{t , as discussed in more detail later. 

The complete hard contribution 

The complete G function, including soft pre-evolution and vertex contributions, 
is given as 

G hard (*+,x-, S ,&) = FiL(x + ) F^x-) (41) 

: f 'dz + dz-EU(z + ) E^{z-) a{Uz+z-x+x-s,Q 2 ) 
Jo 



jk 



1 / b 



■ cxp 



h2 



An\ NN (l/(z+z-)) "V ±Xnn(1/(z + z-)) 

with \nn(0 = + a' soft \n£. The variables x + , x~ refer to the soft nucleon con- 
stituents participating in an elementary interaction. The variables z + , z~ refer to 
the first hard partons, the ladder ends. The function E^ oft has two contributions, 

< ft (z) = E* om (z)+ E* om (z) . (42) 

The first term, representing a soft "pre-evolution", is given as lmT so{t (s /z,t = 0), 
up to some coupling modifications, explicitely 

K m (*) = * _a "* £soft ( o) (*) - K m (*) = ^ ^oft(o) (*) > (43) 

with 

^soft(o)(^) = 7soft^ El m (z), E° om (z) = 7 soft (1 -w)(l- zf, (44) 
with 7 soft = 87Ts 7part7. and 

EL m (z)= j\^ s °P(0(l- Z f, (45) 
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with some parameters 7 and (5. 

The second term, representing a direct coupling to a valence quark, and which 
has therefore only a quark contribution, can be written as 

^soft(i) (x,z) = z- a - E" som (zx,x), (46) 

with 

E q (r r) - r 1 x a ^ ^ + a ™ - (47) 

So this function depends also on the momentum fraction x. We neglected the 
small hard scattering slope Rl ard compared to the Pomeron slope A soft . 

We call E soft also the " soft evolution", to indicate that we consider this as sim- 
ply a continuation of the QCD evolution, however, in a region where perturbative 
techniques do not apply any more. E 3 so{t (z) has the meaning of the momentum 
distribution of parton j in the soft piece. More details can be found in [26]. 



C Appendix on Parameterizations of G Functions 

The calculations of partial cross sections in our multiple scattering theory involve 
high dimensional integrations, excluding numerical techniques. Fortunately, the 
numerically determined G function (the sum of soft and hard) can to a high 
precision be expressed as 



b 2 

G(x + ,x~,s,b) = (x + x-)- a ^ Vac, sP * (x + x-f G i (sx + a0 7G * fe2 exp(-— ). (48) 

i=l Ut 

This form is inspired by the fact that G is a sum of two contributions, soft [i — 1) 
and hard [i = 2), with quite different exponents (3' G .. The fit parameters are, how- 
ever, chosen such that the sum of both, soft and hard, is best reproduced. The 
above form provides in fact excellent fits of the x + , x~, and also the b dependence 
of G. Concerning this b dependence, one has 

(sx + x-y a ' b2 exp(-|^) « exp(-^-), (49) 

with 

\ Gl = S Gl + 5 2 GilG M^ + ^), (50) 
indicating that the width of the b dependence increases logarithmically with en- 
ergy. 

Using the analytical form of G, many of the multidimensional integrals can be 
done analytically. For a given s and b, one may write 

2 

G(x + 1 x- 1 s 1 b) = ^a l (x + x-y% (51) 
i=i 
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with s and b dependent coefficients a*, We also use the notation f3 s = /?i and 
/3# = /3 2 . Taking into account screening via the exponents e s = s 1 and e H = e 2 (see 
section H), we have 

2 

G(x + , x~, s,b) = J2 a i (x + f^ (Zp) (x-y^ (ZT \ (52) 

i=l 

The first task is the calculation of the function $, representing summations 
over all uncut contributions corresponding to a given cut configuration, which we 
need to know to be able to calculate partial cross section. For nucleon-nucleon, 
one obtains 



oo oo 



^ NN ^ + ,x-,s,b) = ( X ^ x ~)^j2T,^^9m+jP2)gm+jP2), (53) 

with a< = -a i r( / 9 i )r( / 9' i )(x + )A(x-)^, and # = fa + e^Zp) + 1, # = A + ^(^t) + 1. The 
function g is given as 

/ \ r(l "I - Q^remn) 

g(x) — — 

r(l + a remn + x) ' 

$ NN is almost an exponential function. So to regularize the theory, we replace 
g(0i + j/3 2 ) by g^PiYg^Y , to get an exponential-type function. It is not precisely 
the same, but we take this freedom, because even small deviations from an 
exponential form of $ may lead to drastic consequences (big oscillations), not 
at all being physical. All these considerations can be generalized to nucleus- 
nucleus (AB) collisions, where we have to evaluate §ab- One gets 

AB 

$ AB (x proj ,a; targ , S ,6) = J]exp(-G fc « (fc) ,x; (fe) , S ,& fe )) (54) 

k=l 

A B 

X n^ + ) Qr ° mn @ ( 1 _ X i) T[( X l) aiemn ®( X j) ©I 1 " Zj), 
i=l j=l 

with 

2 



G k (x + ,x~ } s,b) = J2 ®i ( x+ f' i x ~f' , (55) 



with 



i=i 



k r(# )r(i + « remn ) r(/3'*)r(i + «re mn ) 



a,- = a,- r-, loo) 

r(l + a remn + /f) r(l + a remn + /3'") 

/Qjf = ft + ei{Z P ) + \, (57) 

= 13} + £i(Z r ) + l, (58) 



with 



af = a Gi exp(-^-)/ G i+^, (59) 
A" = flk + Ta^-Op-t- (60) 
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More details can be found in [26] . 



D Appendix on Markov Chain Techniques 

As discussed earlier, the function Q (the integrand of partial cross sections) can 
be interpreted as probability distribution, which allows to apply the Monte Carlo 
technique. So we want to generate configurations {m, X + , X~} according to 

n^\m, X+, x ) = W\— ] W G k (x^, x^, s, b k ) <$> AB (X^°\X^, s, b) , (61) 
k=i I k ' fi=i J 

with 

b = {b k }, m = {m k }, X+ = {xij , X- = {x~ k J , X*«* = {af*} , X** = {xf^} , 

(62) 

Here, b is the impact parameter between the two nuclei, and b k is the transverse 
distance between the nucleons of k th pair, m k is the number of elementary in- 
elastic interaction for the nucleon-nucleon pair k, x k and x k are the light cone 
momenta of the /i th interaction of the pair k. The arguments of $ are the number 
of inelastic interactions m and the momentum fractions of projectile and target 
remnants, 

where n(k) and r(k) point to the remnants linked to the k th interaction. In the 
following, we perform the analysis for given s and b = (bo, h, so we do not 
write these variables explicitly. Furthermore, we suppress the index AB. For any 
given configuration the function Q can be easily calculated, using the techniques 
developed earlier. 

Since Q(m,X + ,X ~) is a high-dimensional and nontrivial probability distribu- 
tion, the only way to proceed amounts to employing dynamical Monte Carlo 
methods, well known in statistical and solid state physics. We first need to 
choose the appropriate framework for our analysis. So we translate our prob- 
lem into the language of spin systems [38]: we number all nucleon pairs as 1, 
2, AB and for each nucleon pair k the possible elementary interactions as 
1,2, m k . Let m max be the maximum number of elementary interactions per 
nucleon pair one may imagine. We now consider a two dimensional lattice with 
AB lines and m max columns, see fig. I28L Lattice sites are occupied (= 1) or empty 
(= 0), representing an elementary interaction (1) or the case of no interaction (0), 
for the k th pair. In order to represent m k elementary interactions for the pair 
k, we need m k occupied cells (l's) in the k th line. A line containing only empty 
cells (O's) represents a pair without interaction. Any possible interaction may be 
represented by this "interaction matrix" M with elements 

m kn e{0,l}. (64) 
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Figure 28: The interaction lattice. 

Such an "interaction configuration" is exactly equivalent to a spin configuration 
of the Ising model. Unfortunately the situation is somewhat more complicated 
in case of nuclear collisions: we need to consider the energy available for each 
elementary interaction, represented via the momentum fractions x\ n and x^ n . 
So we have a "generalized" matrix K = (M, X + ,X ~), representing an interaction 
configuration. Since there are 

AB | 

c = jj rrw (65) 

^ m fe !(m max - m k )\ 

configurations (M, X + ,X ) representing the same configuration (m, X + ,X ), the 
weight for the former is c 1 times the weight for the latter, so we obtain the 
following probability distribution for K = (M, X + , X ~): 

AB ( / x| m max ^ 

m) = n (m T"r n <w * fc - n , s , $ $> AB x**, s , $ . m) 

k=l l max ' n=l J 

In order to generate according to the given distribution Q (K), defined earlier, 
we construct a chain of configurations such that the final configurations 
fT ( ' max) are distributed according to the probability distribution Q(K), if possi- 
ble for a t max not too large! Let us discuss how to obtain a new configuration 
K(t+i) — £ f rom a given configuration = K. We use Metropolis' Ansatz for the 
transition probability 

p(K, L) = prob = L \ K {t) = K) (67) 

as a product of a proposition matrix w(K,L) and an acceptance matrix u(K,L), 
where we use 

in order to assure detailed balance. We use 

«,(*,£)=( J . ^'• L ^ 1 , (69) 

v ' otherwise 
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where d(K, L) is the number of lattice sites being different in L compared to K, 
and where Q is defined by the same formulas as fl with one exception : <3> is 
replaced by Ut=x IT=T <W (l - 4J i 1 ~ *J • So we get 



AB 



fc=i 

II Sm ^ G ( X km x km s > h) (1 - x+ ) (1 - x^ n ) }. 



fc=l 

x 

n=l 



The above definition of w(K, L) may be realized by the following algorithm: 
d choose randomly a lattice site (fc, n), 

d propose a new matrix element (m kn , x\ n , x^ n ) according to the probability 

distribution p(m kn , x\ n , x^ n ), 

where we are going to derive the form of p in the following. From eq. 1701 . we 
know that p should be of the form 

, + -\ ,f G(x + ,x~,s,b) (1 -x + ) (1 - x~) if m = 1 
p(m,x + ,x )~m W ^ if m = ' ( } 

where m = m max — m is the number of zeros in the row k. Let us define m as the 
number of zeros (empty cells) in the row k not counting the current site (fc,p). 
Then the factor m \ is given as m ! in case ofm/0 and as moK^o + 1) in case of 
m = 0, and we obtain 

p(m, x + , x') ~ (m + l)5 m o + a;", s, b) (l - a; + ) (l - x~) 5 mX . (72) 

Properly normalized, we obtain 

p{m,x^,x ) = p d m0 + (1 -Po) <W, (73) 

where the probability p of proposing no interaction is given as 

Po = _ — r— ; 7 r , (74) 
m + 1 + 

with % being obtained by integrating G(l — x + )(l — x~) over x + and x~, 

x(s, 6) = / dx + dx~G(x + , x~, s, b) (l — x + ) (l — aT") . (75) 

Having proposed a new configuration L, which amounts to generating the values 
n^kn, xt n , Xk~ n for a randomly chosen lattice site as described above, we accept this 
proposal with the probability 

u(K, L) = min (ziz 2 , 1) , (76) 
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with 

z ~ n(L) z - W{L ' K) (77) 

Since K and L differ in at most one lattice site, say (k,n), we do not need to 
evaluate the full formula for the distribution Q to calculate z x , we rather calculate 

_ n*»(L) 

with 

Q kn (K) = {(m + l)5 mkn0 + 5 mknl G k (xt n ,x- kn ,s,b k )) 

x expf- £ G^lx^ln)} (79) 

\ £ linked to k / 

x (^)°—e(^)e(i - (^)°— e(x^)e(i - x%*), 

which is technically quite easy. Our final task is the calculation of the asym- 
metry z 2 . In many applications of the Markov chain method one uses symmetric 
proposal matrices, in which case this factor is simply one. This is not the case 
here. We have 

with 

n k n (K) = p(m kn ,xt n ,x k : n ), (81) 

which is also easily calculated. So we accept the proposal L with the proba- 
bility min(ziz 2 , 1), in which case we have K^ t+V) = L, otherwise we keep the old 
configuration K, which means K^ 1 ^ = K. 
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